function theta_star_hat = BS_CI(bs)


% Applying Subsampling method from Politis Romano and Wolf's Chapter 3 
% for stationary time series model

% using block size b = T^.9 = 187. 
% so there are a total of T-b+1 =150 bootstrap samples. 
% so bs =1 : 150

%estimates_star is the BS estimate. 

thisFile = mfilename('fullpath');
thisDir = fileparts(thisFile);
cd(fileparts(fileparts(thisDir)));

% Load the data from data/processed/OilData.mat
data0_true = open('data/processed/OilData.mat');
dataQ_true = table2array(data0_true.OilData.Production)*30/1000;     % dataQ is in million barrels per month
dataP_true = data0_true.OilData.Price;
block_size = floor(length(dataP_true).^0.9); 
I = size(dataQ_true,2); 

%load the theta_MLE and use it for initial value in MLE estimation 
% from ~/output/estimates/ESTIMATES.MAT (this file is provided in this
% replication folder)
load('output/estimates/ESTIMATES.mat'); % loads as theta_mle

% from now on, theta_mle is the initial value
	

% step 0. call the bootstrap sample      
     dataQ = dataQ_true(bs:bs+block_size-1,:); 
     dataP = dataP_true(bs:bs+block_size-1,:);
     T     = size(dataQ,1);
              
% step 1. NLLS estimate of trend parameters

tempfun = @(x,tdata) -exp(repmat(x(1:I),T,1)-x(2*I+1)*tdata)...
    + exp(repmat(x(I+1:2*I),T,1));
tdata = repmat((1:T)',1,I);
tempfun1 = @(x) sum(sum((dataQ - tempfun(x,tdata)).^2));
options = optimoptions('ga','MaxGenerations',250*(2*I+1));
tau_est = ga(tempfun1,2*I+1,[],[],[],[],zeros(1,2*I+1),[9^999*ones(1,2*I) 1],[],[],options);

%step 2. de-trending 
dataQ = dataQ + exp(repmat(tau_est(1:I),T,1)-tau_est(2*I+1)*tdata) ;
dataQ_means = mean(dataQ);
dataQ_std = std(dataQ);
table0 = [dataQ_means' dataQ_std'];

%step 3: kmeans
%[groupid, ~] = kmeans(table0,6,'MaxIter',10000,'Start',table0(start_obs,:));
ng =6; 
% new classification that combines K-means and Geography
groupid = [4;... %Angola
          1;... %UAE
          2;... %Canada
          2;... %China
          4;... %Algeria
          3;... %Ecudaor
          4;... %Egypt
          2;... %UK
          1;... %Iran
          1;... %Iraq
          1;... %Kuwait
          4;... %Libya
          3;... %Mexico
          4;... %Nigeria
          2;... %Norway
          1;... %Qatar
          5;... %Russia
          5;... %Saudi
          6;...  %USA
          3];  % Venezuela


%step4. estimation 
data1 = [dataP dataQ];
conv_ma = (repmat(groupid,1,ng)==repmat(1:ng,I,1));
tempres = estimate_fun(data1,conv_ma,I,ng,theta_mle);
param_est_new = tempres(1:end-1);

theta_star_hat = [param_est_new(1), ... % beta
                  param_est_new(2), ... % lambda
                  param_est_new(3), ... % muU
                  param_est_new(4), ... % sig2U
                  param_est_new(5:ng+4),... %alv
                  param_est_new(ng+5:2*ng+4), ... % betv 
                  param_est_new(2*ng+5:2*ng+6),... % (a1_tilde,  a2_tilde)
                  param_est_new(2*ng+7), ... %barw
                  param_est_new(end), ... %ulow
                  ];

filename = ['output/bootstrap/BS_theta_estimates/theta_star_estimates_' num2str(bs) '.mat']; 

save(filename, 'theta_star_hat')



end



